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Abstract. 

A global solution of the Schrodinger equation for explicitly time-dependent 
Hamiltonians is derived by integrating the non-linear differential equation 
associated with the time-dependent wave operator. A fast iterative solution 
method is proposed in which, however, numerous integrals over time have to be 
evaluated. This internal work is done using a numerical integrator based on Fast 
Fourier Transforms (FFT). The case of a transition between two potential wells of 
a model molecule driven by intense laser pulses is used as an illustrative example. 
This application reveals some interesting features of the integration technique. 
Each iteration provides a global approximate solution on grid points regularly 
distributed over the full time propagation interval. Inside the convergence radius, 
the complete integration is competitive with standard algorithms, especially when 
high accuracy is required. 


PACS numbers: 31.15.p, 02.70.-c, 02.30.Tb, 33.80.-b 
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1. Introduction 

The numerical solution of the Schrodinger equation ihd'&fdt = HA> plays a key role 
in the understanding of molecular dynamics processes. Molecular inelastic pQ and 
reactive collisions [2 have long been treated using quantal wave packet methods. 
In procedures such as the semiclassical treatment of multiphoton processes [3j 
explicitly time-dependent Hamiltonians have to be used to describe the molecule-field 
interaction. A popular and efficient wave-packet propagation algorithm is the Multi- 
Configuration Time-Dependent Hartree method, which is able to handle a large variety 
of high-dimensional problems [1] . Here, our aim is not to deal with many dimensions 
but to develop an efficient way of treating problems with rapidly oscillating time 
dependencies in the Hamiltonian operator or in the wave functions. Several algorithms 
which are used for tinre-resolved experiments with time-independent Hamiltonians 
can be adapted to the time-dependent case and are described in [5] and references 
therein. However for such methods the length of the integration steps must be reduced 
to handle the high frequency terms. The propagation scheme is then based on the 
decomposition of the evolution operator using small time increments. 

N -1 

U(t, 0) = U ((n + l)Af, nAt) (1) 

n —0 

where A t = t/N and 

U{t + At, t) ~ exp [—(i/h)H(t + At/2) At] (2) 

These step by step integration schemes will be referred to in later sections as 
’’continuous methods”. Second order differencing schemes (SOD) [8], split operator 
methods 7) and the short iterative Lanczos propagation [5] all give cumulative 
propagation errors which are proportional to (At) 3 in ([1]) and also errors due to the 
approximate calculation of the action of exp[— (i/h)H At] on the wave function. Other 
high-accuracy integrators which can be used for multi-dimensional systems include 
the generalization of the sympletic partitioned Runge-Kutta method developed by 
Sanz-Serna and Portillo [9] and the method proposed by Kormann [10] in which 
the Hamiltonian H is replaced by a suitable truncation of the Magnus series. 
Unfortunately, these continuous integrators suffer from limitations which could prevent 
their use in some cases. SOD cannot handle non-Hermitian Hamiltonians, while the 
split operator scheme requires that the kinetic operator does not mix coordinates 
and their associated momenta. The calculation of the Magnus development is only 
tractable if the couplings have separate time and coordinate dependencies m- 

A recently proposed global treatment, the Constrained Adiabatic Trajectory 
Method (CATM), QH H2 Q3J f!4] introduces the Floquet Hamiltonian Hp = H — 
ihd/dt defined in an extended Hilbert space including the time coordinate. The 
dynamical Schrodinger equation is transformed into an equivalent partial eigenvalue 
problem in the extended space and the wave function is forced to conform to consistent 
boundary conditions over a short artificial time extension by using a time-dependent 
absorbing potential. In [14] . the CATM was formulated as a global integrator of the 
Schrodinger equation which does not belong to the category of methods described by 
m and (]2|) and it was concluded that the CATM is well suited to the description of 
systems driven by Hamiltonians with explicit and complicated time variations. The 
method does not suffer from cumulative errors (because of its global character) and 
the only error sources are the non-completeness of the finite molecular and temporal 
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basis sets used and the imperfection of the time-dependent absorbing potential. The 
dynamics of a quantum system can also be obtained in a global way (as in the 
CATM) by using the time-dependent wave operator concept (TDWO) [151 . In this 
formalism the dynamics is formally separated into a simple evolution within a given 
active subspace driven by an effective low-dimensional Hamiltonian plus a secondary 
evolution from the subspace to the complete space. The TDWO concept has also been 
generalized |16l to treat almost adiabatic quantum dynamics. This generalization is 
based on a time-dependent adiabatic deformation of the active space and is useful 
for dynamics which do not escape too far from an adiabatic subspace. However the 
effectiveness of both TDWO and CATM hinges on the capacity to integrate some 
nonlinear differential equations and this is the central problem treated in this work. 
A solution of this problem is essential to preserve the principal advantage of these 
formalisms, i.e. their capacity to give a global solution over the whole interaction 
time. To reach this goal an iterative solution seems appropriate. Such a solution was 
proposed in ref. nzi. using a time-dependent version of the Recursive Distorted Wave 
Approximation (RDWA) [15] together with non-linear transformations such as Pade 
approximants [18]. Nevertheless, this iterative solution is not fully satisfactory and 
often fails because of a small radius of convergence and a high sensitivity to the choice 
of the initial guessed solution. 

The present paper proposes a new iterative solution of the time-dependent wave 
operator equations. This is equivalent to solving the time-dependent Schrodinger 
equation without using the approximation of JlJ and ([2]) at any stage of the calculation. 
In section [2] the basic wave operator equations are presented together with the 
iterative integration procedure. In practice we limit ourselves to the non-degenerate 
case (one-dimensional active subspace) which describes the dynamics issuing from a 
given quantum state. This is equivalent to solving the Schrodinger equation for the 
wavefunction. The iteration procedure involves calculating numerous integrals over 
time of matrix and vector elements. Section [3] explores the possibility of using Fast 
Fourier Transforms (FFT) to compute these integrals. In section [4] we illustrate the 
algorithm by studying the transition of a vibrational wave packet between the two wells 
of a model potential energy surface under the influence of two laser fields. Various 
numerical features are analysed in details. Section [5] gives a discussion and conclusion. 

2. Iterative calculation of the time dependent wave operator 

2.1. Wave operator equations 

Let H be the Hilbert space associated with a molecular system and let S a be a model 
subspace of rank m. It is assumed that the initial molecular state is included within 
S 0 . The orthogonal projector corresponding to the model space is called P Q , with 
Pi = P 0 , Pt = p o , tr(P 0 ) = to. During the evolution, the model space is continuously 
transformed into successive active spaces S(t) whose projectors P(t) are solutions of 
the Schrodinger-von Neumann equation: 

ih^>=[H(t),P(t)]. (3) 

The wave operator associated with the two subspaces S a and S(t) is defined as [T5] 
n(t) = P(t)(P 0 P(t)P 0 )- 1 = U(t, 0; H)(P 0 U(t, 0; PjPo)" 1 (4) 



A short iterative scheme within the wave operator formalism 


4 


where U represents the quantum evolution operator associated with the Hamiltonian 
H(t), i.e. ih dU< ' t Qt' H ' > = H(t)U(t, 0; H). In (J4J), (P 0 PP 0 ) _1 is the inverse of P 
within S 0 - It exists only if the Fubini-Study distance between P a and P is small: 
distps(P 0l P) < \ [IB] , i.e. if the dynamics does not escape too far from the initial 
subspace. In this framework the time evolution can be written as 

/ U(t,0-,H)P o = n(t)U(t,0-,H eff ) r ,x 

l H eff (t) = P 0 H(t)Sl(t) W 

In equation [5] the evolution issuing from the initial state is separated into two terms. 
The first one, U(t, 0; H e ff), describes the dynamics within the model space S 0 . The 
second one can be written as f l(t) = P a + X (£) and it induces transitions from the 
model space to the full space, the off diagonal part (A' = Q 0 QP 0 ) inducing transitions 
to the complementary space exclusively ( Q a being the projector associated with the 
complementary space). It is evident that the factorisation of U(H e ff), which includes 
all the fast variations within S a , is an advantage for the dynamical integration process 
if the model space is well chosen. It can be proven that the reduced wave operator 
X (t) = Q 0 flP 0 satisfies [T51 

= Qo(1 “ *(*))#(*)(! + x (t)) p o- (6) 

Solving this equation is the central point of the wave operator formalism. The solution 
gives f l(t) as well as H e ff(t) (cf equation [5]), from which a partial evolution operator 
can be obtained. 

A continuous integration of this non-linear differential equation has been proposed 
in section VI of the review m ■ Unfortunately, many approximations were necessary 
to make the calculation tractable and the final result was not satisfactory in the 
strong coupling regime. To derive a new global solution, we can use the fact that 
Qo{ 1 - X)t|(1 + A)P 0 = Qoj^XPo to rewrite © as 

Qo( 1 - X(t))H F (t)( 1 + X(t))P 0 = 0 (7) 

where Hp(t) is the Floquet Hamiltonian, 

Hp(t)=H(t)-ih 4 - ( 8 ) 

at 

We now derive an iterative treatment, assuming that eq © is not perfectly satisfied 
at a finite iteration order n, 

Qo( 1 - X^ n \t))H F (t)( 1 + X^ n \t))P 0 = A<">(*) ^ 0. (9) 

We are looking for the increment SX which exactly solves the problem: 

Qo(l - X (n) (t) - SX (n) (t))Hp(t)( 1 + X (n) (t) + SX (n) (t)) = 0. (10) 
Expanding dll gives four terms of different significance, 

- Q 0 SX^(t)ffp(t)SX^(t)P 0 
-Q 0 5X^\t)Hp{t)(l+X^\t))P 0 

+Q 0 (l-X^\t))Hp(t)SX^(t)P 0 =0 (11) 

In (fill) , the first term is the error at the previous iteration defined in 0. The 
second term is quadratic in the increment (5X^"^(t)) and can be neglected; this is 
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an approximation commonly used in such iterative treatments. The third term can 
be rewritten as 


Q 0 5X^H F (l+X^)P t 
= QoSX M 

= Q 0 5XWP 0 H(P 0 + Q 0 X^P 0 ) 

r (n) 
l eff 


P 0 H(P 0 + QoX^Po) - Poih—(P 0 


QoX^Po) 


= Q 0 SX^\t)H^\{t)P 0 , 


with 


=P 0 H(t)(P 0 + X^ n \t)). 


The last term in m needs more detailed attention: 

Q 0 (l - X( n \t))H F 6X^(t)P 0 

= Q 0 H 8X™P 0 - QoX^H SX^ih^-Q 0 X^P 0 

at 


( 12 ) 

(13) 

(14) 


We cannot handle differential equations in which (. H5X (”)) or (X^H) couple the 
unknown components of SX 1 -' 1 ') to all the complementary space spanned by Q 0 . We 
shall thus assume that the non-diagonal elements of H in the complementary space, 
given by Q 0 (H — Hdiag)Qoi are small and make two last approximations, 

Q 0 H(t)Q 0 » Q 0 [H(t )] 

diag Q° 


Q 0 X^ n \t)H(t)Q 0 » Q 0 \x^ n \t)H(t) 


- diag 


Qo 


(15) 


This simplification is important but does not have important consequences if the 
procedure still converges. The neglected terms are gradually taken into account during 
the iterative process and finally equation © is satisfied. In summary, the increment 
defined in (1101) can be computed by solving 

ih?-8X^\t) = AW(i) - SX^mifM) + H%Ut) SX^(t) (16) 


dt 


with H e ff defined in Cl and 


HZ g (t) = Qo H(t)-XW(t)H(t) 


- diag 


Qo 


(17) 


In the following we limit ourselves to the case of a one-dimensional model subspace, 
with P a = |z)(z| and Q 0 = Ylk^i !&)(&!> I*) being the initial state of the system. X ^ 
and SX^ becomes vectors and H e ff(t) is a pure scalar function. Equation 
becomes an ordinary differential equation whose rigorous solution for the component 
(fc|<LY( n )|*) = SX { k n) can be expressed using exponential evolution operators, 


8X[ n) ( t ) = exp 





(«) , 


{t') exp 


I f ( 

ihj o v 


H 


(n) 

diag,kk 




This solution is consistent with the initial value of the wave operator 
fi(i = 0) = P 0 , 

hence -Y( n )(0) = 0. Equations (THU) . (1171) and (fT51) are the central equations of this 
paper which are to be solved numerically. 


(19) 
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Similar equations to those given above were derived in m by expanding the evolution 
operator using successive interaction pictures. However they were not applied 
numerically because of the difficulty of calculating the many time integrals which are 
similar to those present in (1181) . Only adiabatic limits could be calculated, with very 
poor results. Here we show that it is possible to design an efficient implementation by 
observing that the global approach represented by equations © and m is consistent 
with the use of a discrete Fourier basis set to span the whole propagation time interval, 
together with the associated periodic quadrature. These tools give effective ways to 
express the matrices H(t)X^(t) and dX^ n \t)/dt present in A^(t) (equation ©)) 
and to calculate the integrals in dl8l) , as will be explained in section [3] However 
the accurate use of Fast Fourier Transforms (FFT) is only possible if all the time- 
dependent functions are made perodic and smooth. 

The Hamiltonian H(t) of a molecule submitted to a laser pulse is the sum of 
a time-independent part H a which represents the unperturbed molecule and of a 
perturbation V{t) with F(0) = V(T a ) = 0, where T a is the laser pulse duration. 
In equation © this corresponds to X^ n \t = 0) = 0 and thus A ^(t = 0) = 0. 
However at the end of the pulse duration, X ^ (t = T a ) can differ considerably from 
its initial value. To impose the continuity and the initial conditions of X^ n \t) one 
should use the CATM procedure mm- To this end, a time-dependent absorbing 
potential is introduced over an artificial time extension [T a ,T\. In the present case, 
this imaginary potential V abs {t) destroys all the components of the wave operator (and 
the wavefunction), except the one corresponding to the initial state (i|f2(i)), 

V abs (t) =-iV opt {t)J2 |fc><*l (20) 

k^i 

where V op t(t) is a real positive function localised on the interval [T 0 ,T]. This 
supplementary term produces the boundary value X^ n \T) —> 0 and hence the 
periodicity of SX^ and (equation ©): 

6XW(t = 0) = dX (n) (f = T) = 0, 

A00(f = 0) =AW(f = T)=0. (21) 


2.3. The problem of nested integrals over time 


The n th order iteration requires solving the equations ©, El and El and 
the passage to the (n + T) th order, via the incrementation of X using X^ n+l ' 1 = 
+ SX an d the construction of A( n+1 ) (cf equation ©). We want to 
use the same discretization at regularly distributed grid points to describe all 
the time dependencies. Suppose that the molecular basis set is composed of N v 
orthonormal states with P a = |*)(*| and Q a = Y^k^i |&)(&| an d that the discrete 
time-grid which spans the time-interval [0,T] is composed of N discrete time values 
tj = jT/N , j = 0, ...,N — 1. The calculation of SX^ n \t) in El calls for 


the calculation of N v — 1 integrals / Q tj (^H^ g 


- H. 


(n) 


) dtf, 


where k takes all the 


N v — 1 index values in the complementary space, followed by N v — 1 integrals 


f Lj A 
Jo 


(n). 


1 (t')e rl eff) aL fa/ xhis is obviously a significant task, which 

must be undertaken at each interation order ( n ). The calculation of each integral 
must be fast and very accurate. The selected procedure is explained in the next 


h. 


(n) 
diag,kk 


-H' 


<»)' 

eff. 


dt 
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section with a preliminary discussion about integrating an arbitrary complex function 
of time. 


3. Discrete Fourier integration at the n th iteration order 

In this section we explain how to compute efficiently a numerical definite integral of 
a time-dependent function described on regularly distributed sampling points. The 
function is not necessarily continuous at the boundaries of the grid. This preliminary 
discussion is important in facilitating the numerical implementation of equation (1181) . 
We require that the procedure should determine with the same accuracy, not only the 
integral at the final time value t = T (which is easily obtained using the standard 
Fourier quadrature rule) but also its value for each of the N intermediate discrete 
time values tj. 


3.1. Continuous integrals and numerical FFT-integrators 

In the present subsection the Fourier transform F(v) is defined as 


/ -+-CO 

f(t ) exp(—i2TTi't)dt. 

-OO 

Let I(t) be the integral 

m= r f (t')dt' 

J — OO 

/ +oo 

f(t')h(t — t')dt' 

-OO 

where h(t) is the Heaviside step function, 


( 22 ) 


(23) 


h{t) = 


0 t < 0 
1 t > 0. 


The Fourier transform of h(t) is the distribution 

„ M = 1 (24 ) 

PV denotes the Cauchy principal value which becomes relevant when the above 
expression is integrated. The convolution theorem gives 


I(t) = f*h(t) 

r+oo 

= / F(v)H(v) exp(2i7rz/)d^. 


Using H(v) from equation (l24l) then gives 


f* +oo 


F(u)e 


2i-KlAt 


-dv. 


(25) 


(26) 


In molecule-field interaction problems, the field term is frequently an oscillating 
function of time. Let us consider for a moment the special case of an oscillating 
function of the form f(t) = g(t) exp(i27ri>f). Its Fourier transform is the translated 
function 


F(v) = G{v-v), 
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where G{y) is the transform of g(t). Using this F(v) in (1^1) and using the variable 
change v' = v — v, we obtain 


/■* 1 

I(t) = / g(t) exp(i2niH)dt = —G{—v) 

J — OO " 


_ glZTVVt py 

2n 


G{-v'y 


-dv'. 


(27) 


When the frequency v is very high, we can neglect the Fourier coefficient G(—v) and 
approximate the denominator u — l/ by u, giving the approximation 

i g(t) 


g{t') exp{2i'Kvt')dt' exp[i(27rz/f — 7 t/2)] 


2 TTV 


(28) 


We can now go back to the general case of equation (E^l) to derive a numerical 
FFT-integrator. The detailed derivation is given in |Appendix A| This appendix 
explains how equation (1261) can be implemented for periodic or non-periodic functions 
known as N sample values distributed over the interval [0,T]. The integration 
algorithm needs only two FFT, with a computational cost scaling as 21V log TV to 
obtain the N values of I (tj), with tj = jT/N,j = 0,1,..., N — 1. Here we recall only 
the main result, 


I{tj) = FFT-\p.FFT(f)) + ± x -±=mu 0 (/). 
where fi.FFT(f) is a component by component vector multiplication, with 


T 


pt = 


J.Z o — i _ i 

27T t c — -L, . . . 2 1 

i T £ _ N 


2tt t-N 


N 
2 


(29) 


(30) 


In addition, the t = 0 component of (p.FFT(f)) in (1291) must be replaced by the 
following quantity: 


JV-l 

a = - E V- FFT t (/)■ 

e=i 


(31) 


In the next subsection we perform some numerical tests using the algorithm derived 
in appendix |Appendix A| with rapidly oscillating functions. 


3.2. Numerical test of the FFT-based integration algorithm 

We study the arbitrary case of a function composed of six gaussians multiplied by 
complex exponential functions. The gaussians have various widths and are distributed 
on the interval [0, T = 180] to give the function 
6 

m = £ exp [—a,j(t — tj) 2 ] enp(i2irXjt/T). (32) 

3 =i 

The various parameters appearing in (1321) are given in table [T] The aj are chosen so 
that the modulus of the function f(t) is equal to or smaller than 10~ 20 at the two 
boundaries (see figure [1]). In this example, the FFT procedure uses the function 
f(t) defined by (IA.6I) . This modified function is characterized by a spectrum with a 
central component £ = 0 equal to zero. This can be seen in figure [2j which shows the 
modulus of the spectrum of f(t). As the functions converge rapidly to zero at t = 0 
and at t = T, the continuity of the function f(t) (assumed to be a periodic function 
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Gaussian 
index j 

dj 

tj 

Xj 

1 

6.790 

27.000 

12.000 

2 

3.819 

36.000 

135.600 

3 

1.018 

90.000 

1.750 

4 

1.591 

108.000 

154.700 

5 

2.118 

135.000 

3.250 

6 

3.310 

144.000 

18.150 


Table 1. Parameters used in equation (f32t . 



time 


Figure 1. |5Re/(£)|, see J32D and table [I] 



Figure 2. Absolute value of the spectrum of /(£) corresponding to dA.6l> and 
(l32t . The function is given as a sample of N = 2048 values equally distributed 
on the interval [0,T]. 
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by the FFT procedure) and also the continuity of the first two derivatives are ensured 
without the need to extend the time interval or to use an asymptotic polynomial (eq 
(IA. 101) 1. Figure [2] confirms that the functions are well represented by the Fourier 
expansion; the components omitted in the spectrum are smaller than 10 -13 . 

We first compare the speed of convergence of two algorithms, the numerical FFT- 
integrator (this work) and the standard Simpson rule. As the exact values of the 
integrals Ij = I[tj) are not known, we perform two series of calculations by doubling 
the number N of sampling points, comparing for each of the two integration procedures 
the result for 2 N points with that for N points. The convergence factor CF n is 
calculated by using the N points (with even indices) that the 2N grid has in common 
with the N grid. At each of the N sampling points the absolute value of the difference 
between the two results for the integral is calculated and the convergence factor is 
defined as their maximum value, namely 

CF JV = Sai(|/g^-/f)|). (33) 

j=0 J J 

The corresponding numerical results are shown in figure [3] As the number of sampling 



Figure 3. Convergence factor 11331) for the FFT algorithm (■) and the Simpson 
algorithm (A). 

points increases, the FFT convergence factor decreases rapidly by many orders of 
magnitude up to a critical N value, Nfp T = 1024. Beyond this point the convergence 
factor reaches a constant plateau with a precision close to machine precision (10~ 14 ). 
This is the usual behaviour expected from the Nyquist-Shannon theorem [26j when 
using Fourier basis expansions. The Simpson rule results show a strongly different 
behaviour. The convergence factor decreases regularly up to a minimum value for a 
critical value , which is about 100 times larger than for the FFT case. Moreover 

the minimum convergence factor reached is about two orders of magnitude greater 
than for the FFT results. After this minimum value, the convergence factor increases 
again. In this new regime the roundoff errors become predominant. Since no exact 
result is available, we next compare the results obtained with the best precision for 
each method, to be sure that the two algorithms are consistent and converge to the 
same value. This corresponds to N = 1024 for the FFT algorithm and to about 
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N = 150000 for the Simpson algorithm. Figure [I] presents for each discrete time 
value the modulus of the difference between the two integrals calculated with the 
FFT algorithm and with the Simpson rule, 

\lFFT {tj) _lSim P{t .)\ ( 34 ) 

This figure reveals that the difference on the whole interval [0,T] is never larger than 
10~ 13 . Similar results, not presented here, are obtained for examples in which it 
was necessary to impose a continuity of the two first derivatives by using a 6 th order 
polynomial (see |Appendix A.21. These applications prove that the FFT integrator 



Figure 4. Difference between the FFT and the Simpson result: | I FFT (tj) — 
/ Slmp (tj)| (full line) and the modulus of the FFT result \I FFT (tj)\ (dashed line). 

possesses the expected qualities. Figure 0] shows that the precision of the results 
is very high and virtually independent of the discrete instant tj taken as the upper 
bound. Moreover, figure [3] reveals that the number of sampling points necessary 
to give convergence for the FFT is 100 times less than that for the convergence of 
the Simpson algorithm. Finally the operations necessary for the FFT algorithm are 
reduced to the calculation of two FFT, the first being made in any case during the 
propagation scheme m- For the selected example the CPU time associated with the 
FFT procedure is about 20 times smaller than that associated with the Simpson rule. 
However, the degree of advantage of the FFT approach is closely related to the shape 
of the spectrum of f(t). In the present cases as well as for the Schrodinger equation 
integration procedure treated in the next section, the spectrums are narrow and are 
thus well represented by a Fourier basis set. 

4. An asymmetric double-well transition experiment 

The illustrative example is that of a model diatomic molecule submitted to two laser 
pulses. Throughout this section we use arbitrary units with the convention h = 1. 
The various numerical parameters have been adjusted so that the relative orders of 
magnitude of the vibrational level spacing, the strength of the dipole couplings and 
the pulse duration produce realistic dynamics with significant transition probabilities. 
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Two potential curves S' = 1,2 are defined, which can refer to the two first electronic 
states of a one-dimensional vibrational Hamiltonian in the framework of the Born- 
Oppenheimer approximation. These two curves are shown in figure [5] The lower 
surface is a double well similar to those used in references [281 : 55] , 

ei (R) = -5R 2 + 0.5R 3 + R 4 . (35) 

The upper surface is a single quartic well, 

e 2 (R) = 0.2i? 4 . (36) 

The R coordinate is chosen as a relative position centred on the barrier of e\. All the 
eigenstates are pure bound states and were calculated using a Fourier grid method on 
the radial coordinate [30:, [31] . They are shown in figure [5] These model potentials are 
also similar to those used to describe, for instance, effective isomerization problems 
along a reaction coordinate [35] • The evolution of the vibrational evolution operator 



R 


Figure 5. The two potential energy curves and some of the first unperturbed 
vibrational eigenstates of the first potential surface (full lines) and of the second 
one (dashed lines). The potentials are quartic polynomials defined as: e\(R) = 
-5R 2 + 0.5H 3 + R 4 ; e 2 (R) = 0.2R 4 . 


is driven by the following equation within the framework of the dipole approximation: 


ipJ(L0) = 


IN ' 


£l ^ 


-fh a .E(t) 

£2 


U(t, 0) 


(37) 


Tn is the relative kinetic energy of the two atoms, i.e. T jv = with an 

arbitrary mass equal to 10. £ 5 = 1,2 are the two energy surfaces. There is no rotation 
in the model and the dipole transition moment ft is assumed to be constant along the 
radial axis R and equal to 1. 

The laser field amplitude projected on the ft direction is defined as the sum of 
two pulses with gaussian envelopes 


m = 


3= 1 


cos (uj(t — Tj )) exp — 



(38) 
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Vibrational states (v = 1,5 = 1) and (v = 6, S' = 1) have in common a strong 
overlap with state (v = 7, S = 2) of the excited potential curve. We have therefore 
chosen the carrier frequencies as aq = E(v = 7, S = 2) — E(v = 1, S = 1) ~ 9.98449 
and u >2 = E{v = 7, S = 2) — E(v = 6, S' = 1) ~ 4.77725. This choice induces 
resonant transitions between the first and the sixth eigenstate of the lower surface: 
uji — UJ 2 = E(v = 6, A = 1) — E(v = 1,S = 1). This will ensure strong transitions 
between the two asymmetric wells and a rather large depletion of the initial state. 
Calculations were made using n = 3.90 and T 2 = 4.50 and two different electric field 
amplitudes and centres, 

Ex = 0.05, Ti = 23.5 

E 2 = 0.08, T 2 = 21.5 (39) 

in the first example and 

Ex = 0.09, Tx = 22.5 

E 2 = 0.05, T 2 = 21.5 (40) 

in the second example. The corresponding electric fields are represented in figure [6] 
Those choices for the laser parameters place the system in a non optimal STIRAP 
configuration [33]. In both cases the initial state is the first eigenstate of the first 
surface (v = 1, S = 1). In terms of the wave operator, the first global iteration starts 
with a guess chosen as f ^(t) = |u = 1 ,S = l)(v = 1 ,S = 1|, i.e. X^lft) = 0. 
A vibrational basis set made of the 30 first vibrational eigenvectors of each surface 



Figure 6. The electric field E(t) in the first example of 113911 (left frame) and in 
the second example of (1401) (right frame). 


(i.e. a total molecular basis set of N v = 60 states) is sufficient to obtain converged 
results. This molecular basis set is small but here the principal aim is to test 
the time propagation algorithm. Here we need N t = 4096 Fourier states | p) with 
(t\p) = exp(*27rpt /T) (or equivalently N t = 4096 time grid points) to describe the 
time evolution correctly. In what follows the vibrational states are numbered from 
v = 1 to v = 30 for the first surface and from v = 31 to v = 60 for the second 
surface. The total basis set describing the extended Hilbert space is then composed of 
N = N v x N t = 245760 states. In this framework, the dynamics is calculated by using 
the iterative scheme explained in section [2] (equations dSJ), (fI31) . 1171) and HID with 
f2(”+i) = P a + X( n+1 1 and V(" +1 ) = + SX^ n \) The integrals app earing in (fl8|) 

are calculated using the numerical integrator presented in section [3] and |Appendix A| 
This scheme is used to calculate all the components of the off-diagonal part of the wave 
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operator for each of the N t discrete t k values: Xy n \t k ), v = 1,..., N v , k = 1,..., N t , 
at each iteration order (n). 



Figure 7. Convergence factor 03 with the parameters defined in (1391) (■) and 
those defined in |40I (C\) with respect to the iteration number. 

Our main concern is the speed of convergence of the global iterative scheme. We 
define the following convergence factor: 

Ef=iEflil^" +1) (4)| 2 ’ 

with SX( n ' ) defined in m- Figure [7] shows that the two calculations converge. A 
plateau of about 10~ 17 for the first example and of about 10 -16 for the second one is 
reached at the 22 th iteration. Between the 10 th and the 22 th iteration, the accuracy 
is improved by one order of magnitude at each iteration. The convergence factor of 
equation m adds the residues 8X^ n \t) over the whole propagation interval. This 
means that the plateau seen in figure [7] indicate a very high accuracy at each time 
grid point. It is essential to note that, by contrast with a continuous propagation 
scheme, in which errors are accumulated at each time step, the global scheme does 
not accumulate the errors between two iteration orders. Nevertheless our iterative 
procedure possesses, as does every iterative treatment, a finite radius of convergence. 
The small increase of the convergence factor observed in figure [7] during the first four 
iterations for the second example is the consequence of a large Fubini- Study distance 
between P a and P in this case (the survival probability at the end of the interaction 
|(u = l|'F(t/inaz)}| 2 is 0.2129 in the first example but only 0.02864 in the second one). 
This increases the distance between P a and P(tf i na i) in this last case. This can be seen 
in figure [8] which illustrates the progressive separation between the active subspace S 0 
and the dynamic subspace S(t). In the non-degenerate case, the Fubini-Study distance 
is simply 

distFs(Po,P(t)) = arccos(|| (uil’F(f)) ||) = arccos 

The right-hand side is obtained using the normalization properties of the wave 
operator. The final distance for the second example of (l40l) is clearly larger. The 


l|n(t)ll)' (42) 
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Figure 8. Evolution of the Fubini-Study distance distpsiPo, P(t)) during the 
pulse defined in J39D (solid line) and (1401) (dashed line). The distance has been 
computed using the wave operator at iteration 22. 


physical reason explaining this strong difference is that example pH) is closer to 
an efficient STIRAP situation than example (1551) (but still non-optimal, the final 
transition probability to state \v = 6) being here only 0.22, because of many non- 
adiabatic transitions to other vibrational states). 

The global character of the integration is evident in figures [9] and [TO] which 
present the transition moduli to the different vibrational states at the final time 
|(v|'I'(")(i/maz)}| at four different iteration orders. It is interesting to note that in 
figure [9] an approximate vibrational distribution close to the exact one is obtained 
as early as at the fourth iteration. The small convergence difficulties during the first 
iteration orders in the second example are evident in figure [101 The 4 th order shows 
some transition amplitude larger than 50, very far from the final results. In this 
case more iterations are needed to achieve a stabilized distribution. Probabilities 
can thus exceed one in the global scheme if the results are observed at an early 
intermediate order, before the completion of the iterative process, such as in Fig. 
m There is no fixed norm conservation in the structure of the algorithm but the 
norm conservation is recovered after convergence. Figures [TT] and [12] also illustrate 
the global character of the time-dependent propagation. In these figures we present 
the evolution of the transition amplitude to state \v = 37). This particular state 
is strongly coupled with the initial state |u = 1). Figure [TT] shows that each 
iteration order produces a projected amplitude which converges towards the exact 
solution over the whole interaction time interval. At the fourth iteration, the difference 
\(v = 37|'p( ra=22 )(t) — ^^(9)1 i s already lower than 10 _1 . At the 9 th order, the same 
error becomes lower than ICR 3 . Figure [12] reveals the difficulties due to the large 
Fubini-Study distance between P a and P(tfi na i) (here the final survival probability 
is |(i> = 1\^{tfi n ai)\ 2 = 0.02864). As a consequence, the amplitude of the projected 
wave function onto the state v = 37 at the A th iteration order is about 10 times larger 
than the exact one for t > 30. Surprisingly, the iterative procedure converges again 
between the 4 t?l and the 9 th order. At the 9 th order the error is smaller than 10“ 2 
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Figure 9. Final transition amplitudes to the different vibrational states 
|(u|^( n ) (tfi na i))\ after the pulse defined in H39D , obtained from the wave operator 
at successive iteration orders. + : (n = 2), X : (n = 4), * : (n = 9), • : (n = 22). 



Figure 10. Final transition amplitudes to the different vibrational states 
|(u|\E f ( n ) (tfinal))\ after the pulse defined in (1401) . at successive iteration orders. 
+ : (n = 2), X : (n = 4), * : (n = 9), • : (n = 22). 


everywhere. 

We think that this efficiency is essentially due to the use of a time-dependent 
effective Hamiltonian in the basic equations of the integration scheme m, da 
and HU), especially when it is compared with previous iterative schemes. Figure [13] 
reveals large oscillation of the effective Hamiltonian defined in (f5]> . The imaginary 
part imposes the unitarity of the evolution in the wave operator formulation. A 
contribution to the evolution operator of the form 

/ 1 rt final 

ex P {jf t J $sm(H e ff)dt' 
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Figure 11. The dashed line represents the modulus of the projected wave 
function (v = 3T Id "'! (£)) | as a function of time at the end of the iterative process 
(n = 22) for the first example (parameters given in (1391) 1 ■ The full lines are 
the defects at intermediate iteration orders: |(n = 37|HiC 7l = 22 )(*) — for 

(n = 4) and (n = 9). 



Figure 12. The dashed line represents the modulus of the projected wave 
function |(u = 37|'l/( rl ^(t))| as a function of time at the end of the iterative 
process (n = 22) for the second example (parameters given in (140 h ). The 
thin continuous lines are the defects at intermediate iteration orders : | (v = 
37|^( n=22 )(t) — 4/( n )(£))| for (n = 4) and (n = 9). The deep black line shows the 
survival amplitude |(u = l|4>( n=22 ) (t))\. 
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Figure 13. Evolution of the real part (continuous line) and imaginary part 
(dashed line) of the effective energy, using the field parameters (1101) . 


is present and gives the correct survival amplitude. The magnitude of these oscillations 
(« 4 for the real part) is even much larger than the distance between the first closest 
vibrational eigenvalues of each surface (for example E(v = 2) — E(v = 1) = 1.451). 
Other iterative methods in which these time variations are neglected completely fail to 
converge on this same numerical example. For example, the adiabatic approximation 
would be justified if in equation (1551) the amplitude \{v\H^ ag (t)\v) — (vi\H^(t)\vi)\ 
is large at each instant, |uj) being the initial vibrational state. In such a case the 
introduction of the approximation (051) into equation (051) leads to: 


6X( n \t) 


ihA^ n \t) 




(43) 


Equation (1551) is very similar to the RDWA iterative method which has already been 
used to compute time-dependent wavepacket propagations [T51IT51J . However it cannot 
be used in the present case. It is evident that the strong time-dependency of H e ff(t) 
(see figure 1T51) and the very small survival probability at tfi na i are inconsistent with 
the use of any adiabatic approximation. Actually the use of (1431) in the integration 
scheme produces a very fast divergence of the iterative series caused by the crossings 
of^andF/g}. Another approximate (but more sophisticated) integration scheme 
was proposed in refs HZ! and m- By using a discrete Fourier basis set: | p) 
with (t\p) = exp(*27rpf/T) derived from the discrete time grid basis set by a FFT 
transformation, a generalised RDWA approximation leads to: 


(v',p'\X {n+1) \vi,p = 0) 

{v',p'\ \{H f - #W)XW + h] | Vi ,p = 0) 


{vup = 0\Hly f \vi,p = 0) - (v',p'\H^\v',p') 


(44) 


Introducing (1441) into the iterative scheme also produces a divergence, albeit slower 
than the one produced by (1551) . Undoubtedly this defect comes from the fact that 
equation (1441) only uses an averaged H e f f, represented by the first Brillouin component 
(p = 0\H e ff\p = 0), rather than the exact and strongly time-dependent expression. 
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Finally it is crucial to verify that the global scheme gives the correct results by 
a comparison with some continuous wave-packet propagation algorithms. We first 
compare our results with those obtained for the same system by using a split operator 
[7 procedure by separating the role of the diagonal unperturbed Hamiltonian and 
the role of the laser-molecule coupling. This last contribution is not diagonal and its 
role is approximated by using a Second Order Differencing (SOD) scheme. We also 
calculate the same dynamics using the short iterative Lanczos (SIL) method [8]. To be 
significant, the same vibrational basis set of N v = 60 states is used in the continuous 
propagations (for the SIL propagations, the size of the Lanczos subspace was set up 
to 10). A convergence factor is defined as 



Figure 14. Convergence factors J~c and Tq defined in equations 1451 and dp) 
for the global method (circles), the short iterative Lanczos scheme (triangles) and 
the split-operator algorithm (squares) using the field parameters of equation d39D 
(□, A and o) and of equation J40D (■. ▲ and •), as a function of the number of 
time steps. 


N v 

fc(N t ) = Y I n c ’ Nt (i,t final) - n GA 096 (i,t final )\ 2 . (45) 

i =1 

^G',4096 

is the wave operator calculated at the final time using the global propagator 
with 4096 time grid points after 22 iterations and D c is the wave operator 
reconstructed by using equation after a complete propagation with the continuous 
propapagators with N time steps. To show the convergence of the global method itself 
we have also computed a similar convergence factor between global results, 

JV„ 

F G (N t ) = Y\n G ’ Nt (i,t final )^n G ’ 8192 (i,t final ) 1 2 . (46) 

i= 1 

Figure [IH shows the evolution of Tc and Tq with respect to the number of discrete 
time steps. This figure reveals that the three procedures converge to the same final 
values (same amplitudes and same phases). Figure [TH also confirms that the number 
of time steps needed to converge to a given accuracy is much smaller with the global 
scheme than with the continuous schemes. The convergence is faster for the first pulse 
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than for the second one. The number of time-steps of length dt for the continuous 
integrator is quite large especially if a high accuracy is needed. The CPU time needed 
to compute SX^ (equation (1181) 1 at a given iteration order is of course much larger 
than that required to construct one step in the continuous algorithms. Nevertheless 
the convergence plateau in figure [7] are reached after only 22 iterations. With the field 
parameters of equation (1391) . passing from iteration 22 to iteration 23 indicates that 
the transition amplitudes are converged with 7 stable and significant digits. More 
detailed informations about the CPU time are given in table [2] with varying required 
accuracies. This table shows the ratio between the CPU times required to reach a 
given threshold of convergence using the global or the continuous strategies, using the 
convergence factors Fg and Tc defined in (l45l) and (l46l) . The results of the global 
method have been observed after 22 iterations. The continuous algorithms are faster 


Convergence factor < 

10" b 

nr 8 

io- iU 

10”^ 

io - i4 

Stable digits ~ 

3 

4 

5 

6 

7 

rpCPU /rpCPU 

3.83 

3.11 

0.775 

0.757 

0.211 

rpCPU /rpCPU 

1 G / ^ Split—Op. 

15.3 

6.52 

3.31 

1.62 

0.912 


Table 2. Approximative ratios of the CPU times required to reach a given 
convergence threshold between the global (G) and the continuous methods (SIL 
or split-operator), for example of equation JMJ. The corresponding number of 
significant digits on the largest transition probabilities are also given. 

if a low accuracy is sufficient. The SIL algorithm and the global method then give 
comparable CPU times for intermediate accuracies, while the split operator stays a 
little faster. The global method becomes slightly faster if a high accuracy is desired. 

5. Conclusion 

The theoretical analysis and the simple illustrative examples presented in this paper 
demonstrate that a short global iterative scheme can solve the time-dependent 
Schrodinger equation with an explicitly time-dependent Hamiltonian operator. This 
global method is competitive with a continuous step by step algorithm. The following 
points should be highlighted: 

i) A global integrator does not accumulate discretization errors along the 
time interval and consequently can generate high precision results. The numerical 
experiments have shown that each new iteration order increases the precision by 
one order of magnitude. The structure of the algorithm also allows for a rigourous 
estimation of the error associated with each iteration order by following the quantity 
HpQ — f IHQ, the calculation being stopped when the expected accuracy is reached. 

ii) The use of an FFT procedure to calculate the numerous time integrals appears 
as a very efficient choice which is compatible with producing a highly accurate solution 
to the Schrodinger equation. Nevertheless it is essential that the Fourier spectrum of 
the perturbation should be well reproduced by the Fourier basis set and it is certainly 
more difficult to solve problems with broader spectra. 

iii) The global iterative scheme does not introduce limitations concerning the 
nature or the structure of the Hamiltonian which drives the dynamics. Future 
applications could usefully be aimed at checking how well it can reproduce dissipative 
dynamics. 
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iv) The inability of previous global iterative techniques to determine the time- 
dependent wave-operator for similar problems is partly due to a wrong estimation of 
the effective Hamiltonian acting in the selected model space. In the selected examples 
H e ff(t) exhibits large oscillations which must correctly be taken into account. 

v) In the present applications the starting input for the wave-operator is simply 
X( ra= i) = 0. The iterative nature of the method makes possible a rapid generation 
of new solutions, after a slight modification of the Hamiltonian parameters, by 
simply starting with the wave-operator obtained from the previous solution. Many 
applications are possible, among which we can mention the study of laser control 
processes using exceptional points of parametric Hamiltonians |34i, 05] 06] 07] or of 
vibrational cooling using zero-width resonances [38] . 

The algorithm proposed in this paper appears to be useful for investigating new 
formulations of quantum control in molecular physics or for treating complicated near- 
adiabatic evolutions, for which the use of an effective Hamiltonian can efficiently reduce 
the dimension of the Hilbert space [55] . The global algorithm add the time as a 
new quantum coordinate, which can be a disadvantage to compute the dynamics of 
high dimensional problems over long times. On the other hand the method is not 
dependent on the structure of the Hamiltonian as it can be the case for the split 
operator algorithm or Magnus expansions. The wave operator approach is in principle 
transposable to dynamical processes like dissociation, ionization, dynamics involving 
degenerate eigenvalues or non-adiabatic interactions. For dissipative processes, a 
special difficulty can appear in situations for which the dynamics escapes too far 
from the selected model subspace. In such cases, the algorithm can diverge. The 
best option is to generalise the iterative algorithm of section [2] to handle a degenerate 
active space, by incorporating all the states which are strongly coupled into the model 
space. Study of this generalisation is in progress, as the adaptation of the algorithm 
to the calculation of Floquet eigenstates. 
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Appendix A. Numerical FFT integration 


Appendix A.l. Integration of periodic functions 
Consider equation ESI) . 


m = \m-^ pv 


p+OO 


F(v)e 




-do, 


(A.l) 


where F(y) is the Fourier transform of fit). If the function fit) is assumed to be non¬ 
zero only over a finite interval [0, T], then a discrete finite time-grid can be introduced 
to span this time-interval and also the corresponding finite frequency representation. 
For even N we have: 

jT 

N ’ 


ti = 


j = 0,... ,N — 1 , 
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7 „ iV 

Ua = —, 1=0,...,-1 , 

J rpl j 2 

N 

vn/2 - > 

N 

Vj = -VN-j, j = — + 1, •. ■, A - 1. (A. 2) 

If the time t and frequencies v are selected as one of the points of the above discrete 
grid, a discrete version of the integral can be computed. The continuous Fourier 
transform is approximated by a discrete Fast-Fourier Transform (FFT): 


where 


/ +oo 

f (t)e~ 2i ' KVkt dt = 

-OO 


' —OO 

r T j 1 

/ f{t)e- 2nVkt dt 

I o 


Vn 


FFT k (f) 


N—l 


FFT k (f) = -L Y, f je - i2 * jk ' N 
v A j=0 


(A.3) 


(A.4) 


with fj = f(tj). The integral in equation (l26l) is then apparently approximated by 
the sum 




1 T 
2//V 


FFT k = 0 (/) 


±J_ FFTe(f) c2Mi/N 

2n VN ^ vt 


(A.5) 


However, the passage from the continuous expression (l26l) to the FFT equivalent 
(equation (1A.5[| ) has not been carried out with sufficient rigour, since there is a 
division by 0 for t = 0 in (IA.5|l . We must adjust the term t = 0 before doing 
the inverse transform. The use of a Cauchy principal value and the presence of a 
singularity in (1261) ati/ = 0ifF(j/ = 0)^0 are inconsistent with the use of a Fourier 
series expansion. The equations should thus be modified to take into account this 
discrepancy. However, we can preserve the feature that only two Fourier transforms 
are required to obtain the integral. 

The calculation of the N integrals / (tj). j = 0,..., IV — 1 is made as follows. The 
Fourier transform FFT(f) is first calculated and then a new function f(t) is defined: 

m = f(t)--±=FFT k = 0 (f). (A.6) 


Both / and / have the same discrete spectrum, except for the £ = 0 component, which 
for / is equal to zero, so that the integral of f(t) over the interval [0, T] is equal to zero. 
This is equivalent to using a translation of the ordinate axis. We work temporarily 
with / instead of /. We assume that the spectrum F(v) (equal to zero for v = 0) has 
a Taylor expansion near v = 0: 


F(v) = av + bv 2 + ... 


We are looking for the finite coefficient a = lim Afid, g0 ag apply equation (IA.5I) 

at l = 0 in a consistent way. The divergence produced by the term vn oc \ for £ = 0 
in (IA.5I) has disappeared but the value of a is still unknown. One can guess it by 
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requiring that /(to = 0) = / 0 ° f(t)dt = 0. Equation (IA.5D in this case leads to the 
value 

N-l 

(A.7) 

where p.FFT^f) denotes a component by component vector multiplication, with 

(A.8) 

This value of a replaces the / = 0 component of ( p.FFT{f )) in equation (IA.9I) below. 
So far we have ignored the zero frequency component of the function by working with / 
instead of /. The final result should then be corrected by a linear correction function, 
giving 

I(t 3 ) = FFT~\p.FFT{f)) + ± x J-FFT k = 0 {f). (A.9) 

The first term on the right hand side represents the integral of / at the discrete time 
values tj and the second term restores the correct value of f n 3 f(t)dt. Equations (IA.9I) . 
(IA.8I) and (IA.7D are the basic equations needed to calculate the N integrals I(tj) using- 
two FFT calls (with a computational cost scaling as 2AHogAf). 

The above equations have close similarities with those of the FFT-based 
integration algorithms proposed in ref. [20L 121] in the framework of image resampling. 
The derivation starting from (1231) emphasizes the difficulties introduced by the Cauchy 
principal value term for small frequencies. Here the principal value integral is 
approximated by first cancelling F(v = 0) and guessing its first derivative at v = 0 in 
a consistent and unambiguous way. Our expressions also work well for functions with 
integrals over the interval [0, T] which are not necessarily equal to zero. 

Appendix A.2. Integration of non-periodic functions 

In the algorithm developed in section [2] the smoothness of the wave operator and 
the wavefunction at the time boundary t = T is ensured by the presence of a time- 
dependent absorbing potential (cf. section 12.211 . However the intermediate matrix 
elements of H^(t) and H^ ag kk {t) to be integrated are not necessarily T-periodic. 
FFT-based differentiation and integration algorithms are very accurate but suffer from 
boundary effects (Gibb’s phenomenon). The present algorithm will give accurate 
results only with periodic continuous functions, ideally those with time derivatives 
which are themselves also continuous functions. There are different ways of adapting 
a Fourier-based algorithm to deal with non-periodic signals. For example, extension 
methods have been developed to solve partial differential equations on a complex 
domain [22] [23], to represent surfaces of complex bodies [23] and to treat imaging 
problems [20] . 

Here we use the alternative procedure of Hermite interpolants to ensure the 
continuity of the function /(f) and of the first derivatives [25]. Let /(f) be a complex 
continuous L 2 function defined on a large time interval [0,To], and represented by 
a sample of No equally distributed values /(/?), tj = jTo/N 0 , j = 0, ...,No — 1 
on this interval. The procedure described in (IA.6I) has been applied to /(f) so that 
FFT k =o(f) = 0. As the function /(f) is assumed to be continuous with continuous 
first and second derivatives inside the interval [0, To], the principal problem is to ensure 


H = 


— ) 2?r( 


— - 1=1 


27T t—N 


...£ - 1 
1 = f,...Ar- 1 
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continuity at the boundaries. Let /(0), /'(()), /"(O) and f(T 0 ),f'(T 0 ),f''(T 0 ) be the 
values of the function and of the first two derivatives at the two boundaries. It will 
be impossible to obtain accurate results if these values do not match. One can extend 
the time interval by adding a narrow extra interval [To,T] spanned by an extension of 
the grid tj = jTo/No, j = No, ...,7V— 1 and then constructing on this new interval 
a 6 th order polynomial: 

6 

f(t) = Y,a 3 (t-Toy. (A.10) 

3 =0 

Seven conditions are imposed on the polynomial. The first three are related to the 
continuity of f(t) at the point tj = To and give the first three coefficients a 3 : 

a o = /(T 0 ), 

ai = /'(T 0 ), (A.11) 

0-2 = /"(T 0 )/2. 


The three values on the right hand sides of (IA.Ill) are known exactly if the function 
f(t) has an analytical expression. Alternatively, the derivatives can be approximated 
by a FFT algorithm if fit) is known at a sample of No values. The last four coefficients 
are related to the continuity between t = T and t = 0 of f(t) and to the condition 

fo ftf) dt = 0: 


/(T) = /( 0), 
f(T) = /'(0), 
f"(T) = /"(0), 

Ito f^ dt = °- 


(A.12) 


With A = T—Tq, the last four unknowns are found 


by solving the small linear system: 


A 3 

A 4 

A 5 

A 6 \ 


/ a 3 \ 

3A 2 

4A 3 

5A 4 

6A 5 


a4 

6A 

12A 2 

20A 3 

30A 4 


a 5 | 

A 4 

4 

A 5 

5 

A 0 

6 

) 


\ °6 / 

( /(o) 

- ao 

— a\A — 

a 2 A 2 \ 



/'(0) - ai - 2a a A 
/"(0) - 2a 2 

A /\2 A 3 

\ -a 0 A-ai—-a 2 — ) 


(A.13) 


Note that the last (integral) condition (IA.12I) is not absolutely necessary and can be 
relaxed by first defining the above polynomial and then applying the zero-integral 
condition on the whole interval using equation (IA.6II . 

In summary, the FFT-based integration algorithm is implemented using (1A.7I) - 
(IA.9I) . completed by a preliminary treatment on the additional interval [T 0 ,T] using 

(Iato1) - (ICT1) . 
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